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ABSTRACT 

This  thesis  examined  the  use  of  response  surface  methodology  (RSM)  as  a  parameter 
estimation  technique  in  the  field  of  groundwater  flow  modeling.  Using  RSM,  an  attenq)t 
was  made  to  calibrate  three  hydraulic  parameters  (porosity,  transverse  permeability,  and 
rate  of  recharge)  of  an  existing  two-dimensional,  steady-state  flow  model.  The  model 
simulated  groundwater  flow  for  a  portion  of  landfill  10  located  on  Wright-Patterson  Air 
Force  Base,  Ohio.  The  model  had  previously  been  calibrated  by  grtq)hical  matching 
observed  water-levels  to  predicted  water-levels.  Using  the  parameter  values  from  the 
earlier  calibration  effort  as  a  starting  point,  a  central  composite  design  was  developed  and 
the  simulation  executed  at  each  design  point.  A  residual  sum  of  squares  function  was 
used  as  the  calibration  criteria  and  an  empirical  model  of  the  error  surface  was  developed. 
Of  the  three  hydraulic  parameters,  only  transverse  permeability  had  a  significant  effect  on 
the  response.  The  regression  model  also  indicated  the  response  had  a  high  degree  of 
variability.  A  graph  of  the  regression  equation  revealed  no  local  optima  within  the  design 
region  indicating  the  initial  parameter  estimates  may  not  have  been  warranted. 
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ESTIMATING  GROUNDWATER  FLOW  PARAMETERS 
USING  RESPONSE  SURFACE  METHODOLOGY 


I.  Introduction 


Background 

The  demand  for  municipal,  agricultural,  and  industrial  use  of  groundwater,  and  concerns 
about  groundwater  contamination,  have  steadily  grown  over  the  years.  In  response, 
hydrologists  and  water  resources  personnel  have  searched  for  better  ways  to  understand 
and  manage  groundwater  systems.  This  search  has  led  to  the  development  of  numerous 
groundwater  flow  and  solute  transport  models.  The  majority  of  these  models  are 
simulation  models  implemented  on  computers.  Generally,  these  simulation  models 
include  four  categories  of  mathematical  descriptions: 

1 .  the  relevant  laws  which  govern  the  various  processes  inherent  in  the  system; 

2.  the  hydrogeology  (aquifer  configuration  and  parameters); 

3.  the  external  forces  exerting  on  the  system;  and 

4.  the  initial  and  boundary  conditions  of  the  system. 

While  the  relevant  laws  are  general  principles  that  apply  to  every  groundwater  flow 
system,  the  parameters  of  the  second,  third,  and  fourth  categories  must  be  specified  for 
the  model  to  adequately  characterize  a  particular  system  (Xiang  and  others,  1993:1661). 
Adequately  determining  these  parameters  is  one  of  the  most  difficult  aspects  of 
groundwater  flow  modeling  (Neuman,  1973:1(X)6). 

The  process  of  estimating  the  model  parameters  to  obtain  a  reasonable  match  between 
observed,  site-specific  data  and  model  calculations  is  known  as  model  calibration 
(Walton,  1992:35).  Traditionally,  model  calibration  has  been  done  by  manual  trial-and- 
eiror  and  graphical  matching  techniques.  In  addition  to  being  time  consuming,  these 
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methods  are  often  very  qualitative  and  leave  the  reliability  of  the  model  in  question  (Faust 
and  Mercer,  1980:572).  While  several  automated  techniques  have  been  developed  to 
estimate  parameters  by  optimizing  an  objective  function,  most  practitioners  still  prefer  the 
manual  trial-and-error  method  (Anderson  and  Woessner,  1992:266). 

One  approach  that  may  be  useful  in  performing  model  calibration  is  response  surface 
methodology  (RSM).  RSM  consists  of  a  set  of  techniques  used  in  the  empirical  study  of 
relationships  between  one  or  more  responses  and  a  group  of  input  variables.  These 
techniques  include: 

1 .  Designing  a  set  of  experiments  that  will  yield  adequate  measurements  of  the 
response(s)  of  interest. 

2.  Determining  a  mathematical  model  that  best  fits  the  data  coUected  from  the 
experimental  runs. 

3.  Determining  the  optimal  setting  of  the  experimental  factors  that  produce  tht 
maximum  (or  minimum)  value  of  the  response  (Khuri  and  Cornell,  1987:3). 

Specific  Problem 

The  goal  of  this  research  was  to  determine  if  RSM  could  be  effectively  used  as  a 
calibration  technique  to  estimate  groundwater  flow  parameters. 

Research  Objectives 

To  accomplish  the  stated  goal,  the  following  objectives  were  established: 

1 .  Examine  a  previously  calibrated  model  that  simulates  groundwater  flow  for  a  cross- 
section  of  landfill  10  located  on  Wright-Patterson  Air  Force  Base,  Ohio. 

2.  Ensure  the  simulation  model  implemented  on  available  computer  systems  produces 
results  similar  to  those  obtained  by  the  developer  of  the  model. 

3.  Review  the  available  field  data  associated  with  the  model  and  establish  a  calibration 
criterion. 
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4.  Attempt  to  recalibrate  the  model  using  the  RSM  approach. 

5.  Compare  both  the  effectiveness  and  efficiency  of  RSM  verses  the  manual  trial-and- 
error  method. 

Scope  and  Limitations 

1 .  This  study  used  an  existing  two-dimensional,  steady-state  flow  model  of  a  landfill  to 
demonstrate  how  RSM  could  be  applied  to  the  problem  of  parameter  estimation.  The 
developer  of  the  model  used  a  graphical  matching  technique  to  calibrate  some  of  the 
model's  hydraulic  parameters.  It  was  assumed  the  calibrated  parameter  values  were 
close  to  optimal.  However,  no  quantitative  information  on  the  closeness  of  calculated 
response  to  observed  data  was  available.  It  was  also  assumed  that  the  model  was 
valid  and  adequately  represented  the  physical  characteristics  of  the  landfill  (including 
aquifer  configuration,  external  forces  exerting  on  the  system,  and  boundary 
conditions). 

2.  This  calibration  process  focused  on  three  of  the  most  commonly  adjusted  hydraulic 
parameters  and  used  the  calibrated  parameter  values  determined  by  the  developer  of 
the  model  as  the  starting  point  of  the  investigation. 

3.  The  available  field  data  provided  only  two  points  for  comparing  the  observed  values 
with  the  estimated  response  from  the  model. 

4.  No  information  was  available  on  the  measurement  errors  associated  with  the  observed 
head  values. 
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n.  Parameter  Estimation  and  the  Calibration  Process 


Introduction 

This  chq)ter  reviews  the  parameter  estimation  techniques  used  to  calibrate  groundwater 
flow  models  and  discusses  the  use  of  response  surface  methodology  as  a  possible 
parameter  estimation  procedure.  The  problem  of  parameter  estimation  is  set  in  context  by 
introducing  some  basic  principles  of  groundwater  hydrology,  the  governing  differential 
equations  used  to  model  groundwater  flow,  and  the  process  used  to  develop  groundwater 
models.  The  ch^ter  concludes  with  an  overview  of  RSM  and  how  it  may  be  useful  as  a 
parameter  estimation  technique. 

Groundwater  Hydrology 

Geologic  Framework  and  Terminology.  Any  analysis  of  subsurface  fluid 
flow  must  take  into  account  the  geologic  setting.  This  setting  includes  the  different 
sediment  and  rock  types  in  the  study  area,  spatial  and  temporal  relationships  among  the 
different  subsurface  formations,  and  some  sense  of  the  scale  of  spatial  variability. 
Incomplete  characterization  of  the  geologic  setting  can  lead  to  enors  or  misinterpretations 
in  groundwater  investigations  (Smith  and  Wheatcraft,  1993:6.2). 

Aquifers.  The  aquifer  is  the  most  commonly  studied  subsurface 
formation  in  groundwater  hydrology.  An  aquifer  is  a  geologic  formation  which  contains 
water  and  allows  large  amounts  of  water  to  move  through  it  in  response  to  physical 
forces.  The  porous  media  comprising  the  aquifer  may  consist  of  several  sediment  and 
rock  types.  Regardless  of  its  conqx)sition,  the  portion  of  the  formation  which  is  occupied 
by  solid  matter  is  called  the  solid  matrix.  The  remaining  portion  is  called  the  void  space 
or  pore  space.  This  pore  space  is  occupied  by  water,  or  gases  (mainly  air  and  water 
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vapor),  or  some  combination  of  water  and  gases.  The  zone  of  saturation  is  defined  as  the 
area  of  a  subsurface  formation  where  the  pore  space  is  completely  filled  with  water  (see 
Figure  2.1).  While  the  term  subsurface  water  is  used  to  denote  all  the  water  beneath  the 
surface  of  the  ground,  the  term  groundwater  refers  specifically  to  the  water  in  the  zone  of 
saturation  (Bear  and  Verruijt,  1987:1-4). 
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Figure  2.1  Subsurface  Moisture  Zones  (Bear  and  Vemiyt,  1987:4) 


Most  aquifers  can  be  considered  as  underground  storage  reservoirs  that  receive 
recharge  from  rainfall.  The  rate  at  which  water  flows  into  an  aquifer’s  zone  of  saturation 
is  known  as  rate  of  recharge.  Only  a  fraction  of  the  precipitation  eventually  reaches  the 
zone  of  saturation  and  becomes  part  of  the  aquifer's  groundwater  system.  Figure  2.2 
shows  the  average  annual  hydrogeologic  budget  for  Southwest  Ohio.  While  the  area 
receives  an  average  yearly  rainfall  of  36  inches,  the  rate  of  recharge  is  only  about  six 
inches  per  year. 

Hydraulic  Properties.  To  understand  and  characterize  groundwater 
systems,  hydrologists  have  defined  several  properties  of  a  porous  material.  Ttw  most 
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commonly  used  definitions  are  porosity,  hydraulic  conductivity,  permeability,  specific 
storage,  transmissivity,  and  storage  coefficient. 


Figure  2J2  Annual  Hydrogeologic  Budget  for  Southwest  Ohio 
(OSRI  Report,  1993:F1.7) 


In  the  zone  of  saturation,  the  porosity  is  a  direct  measure  of  the  water  contained  per 
unit  volume,  expressed  as  a  ratio  of  the  volume  of  pore  space  to  the  total  volume. 
Porosity  is  a  dimensionless  number  less  than  one,  although  it  is  frequently  reported  as  a 
percentage.  Table  2. 1  shows  a  range  of  porosities  for  a  several  types  of  geologic  media. 
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Table  2.1  Representative  Values  oi  Porosity  (Smith  and  Wheatcraft,  1993:6.9) 


Sediment  or  Rock  Type 

Porodty 

Sediment  or  Rock  Type 

PoTority 

Clays 

0.40-0.60 

Fractured  igneous  rocks 

0.01-0.10 

Silts 

0.35-0.50 

Basalts 

0.01-0.25 

Fine  sands 

0.20-0.45 

Carbonate  mud 

0.40-0.70 

Coarse  sands 

0.15-0.35 

Dolomite 

0.001-0.15 

Shales  (near-surface) 

0.30-0.50 

Tertiary  limestone 

0.20-0.35 

Shales  (at  depth) 

0.01-0.10 

Paleozoic  limestone 

0.001-0.10 

Sandstones 

0.05-0.35 

Oolitic  limestone 

0.01-0.25 

Bedded  salt 

0.001-0.005 

Karstifred  limestone 

0.05-0.50 

Unfractured  isneous  rocks 

0.0001-0.01 

Chalk 

0.15-0.45 

Hydraulic  conductivity  (JO  is  a  measure  (m/sec)  of  an  aquifer's  ability  to  transmit 
water.  Hydraulic  conductivity  is  a  function  of  both  the  medium  and  the  fluid  and  can 
vary  over  many  orders  of  magnitude.  To  separate  the  effects  of  the  medium  from  those  of 
the  fluid,  the  permeability  (k)  is  defined  as 


where  p  (kg/m  sec)  is  dynamic  viscosity  of  the  fluid,  p  (kg/m^)  is  fluid  density,  and  g 
(m/sec^)  is  the  acceleration  due  to  gravity.  Thus,  permeability  (m^)  is  a  property  of  the 
medium  only  and  describes  how  well  the  medium  transmits  water.  When  permeability 
differs  according  to  the  direction  of  flow,  the  permeability  is  said  to  be  anisotropic. 
Direction  independent  permeability  is  isotropic.  When  permeability  is  anisotropic,  there 
is  always  one  particular  direction,  Xp,  along  which  permeability  has  an  absolute  maximum 
value,  k^.  Somewhere  in  the  plane  normal  to  the  maximum  direction,  there  is  a 
direction,  x„,  in  which  permeability  has  the  absolute  minimum  value,  kg^.  An 
anisotropic  permeability  field  in  two  dimensions  is  completely  described  by  the  extreme 
permeability  values,  k^,,^  and  kg^,  and  the  angle  orienting  the  principal  directions,  Xp  and 
Xg,,  to  the  X  and  y  directions  (see  Figure  2.3)  (Voss,  1984:28).  The  maximum 
permeability  value,  kg^,  is  often  referred  to  as  longitudinal  permeability  while  the 
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minimum  permeability  value,  is  referred  to  as  transverse  permeability.  The 
anisotropic  ratio  of  permeability  is  the  ratio  of  longitudinal  to  transverse  permeability.  In 
layered  geologic  media,  this  ratio  may  exceed  values  as  large  as  100: 1  (Smith  and 
Wheatcraft,  1993:6.9). 


Figure  23  Anisotropic  Permeability  (Voss,  1984:29) 

Specific  storage  (Sg)  is  defmed  to  be  the  volume  of  water  released  from  a  unit  volume 
of  aquifer  per  unit  change  in  hydraulic  head  (see  description  below).  Specific  storage 
(1/m)  is  a  function  of  both  the  compressibility  of  water  and  the  compressibility  of  the 
porous  medium  (Anderson  and  Woessner,  1992:17). 

Transmissivity  (7)  and  storage  coefficient  (5)  are  commonly  used  when  modeling  two- 
dimensional  horizontal  flow  through  an  aquifer  with  a  specified  layer  thickness. 
Transmissivity,  T=  Kb  (mVsec),  is  the  product  of  hydraulic  conductivity  and  the  layer 
thickness;  and  the  storage  coefficient,  S  =  S^b  (dimensionless),  is  the  product  of  specific 
storage  and  the  layer  thickness,  where  b  (m)  is  the  thickness  of  the  layer  being  modeled 
(Smith  and  Wheatcraft,  1993:6.12). 
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Groundwater  Flow. 

Hydraulic  Potential  and  Fluid  Fhix.  Understanding  the  movement  of 
groundwater  requires  a  knowledge  of  the  time  and  space  dependency  of  the  flow,  the 
nature  of  the  porous  medium  and  fluid,  and  the  boundaries  of  the  flow  system. 
Groundwater  flows  through  interconnected  void  spaces,  along  microcracks  between  grain 
boundaries,  and  in  larger-scale  fractures.  Groundwater  moves  in  response  to  differences 
in  fluid  pressure  and  elevation.  The  driving  force  is  measured  in  terms  of  hydraulic  head 
(h),  defined  as 

h  =  -^+z  (2-2) 

Pg 

Here  p  (N/m^)  is  the  pressure  of  a  fluid  with  constant  density  p  (kg/m^),  g  (m/sec^)  is  the 
acceleration  due  to  gravity,  and  z  (m)  is  the  elevation  of  the  measurement  point  above 
some  reference  elevation  or  datum.  Hydraulic  head  (m),  also  referred  to  as  piezometric 
head,  is  equal  to  the  mechanical  energy  per  unit  weight  of  the  fluid.  Groundwater  flows 
from  regions  where  the  hydraulic  head  is  higher  towards  regions  where  it  is  lower. 
Defining  pressure  head  as 


leads  to 


(2-3) 


h  =  h^+z  (2-4) 

where  z  is  the  elevation  head  (see  Figure  2.4).  The  water  table  is  defined  as  the  surface 
on  which  the  pressure  head  is  equal  to  zero.  Contour  m^  of  hydraulic  head  are  used  to 
infer  directions  of  subsurface  fluid  flow  since  flow  will  be  normal  to  the  head  contours  in 
an  isotropic  medium  (Smith  and  Wheatcraft,  1993:6.6). 

Dare's  Law.  The  movement  of  groundwater  is  well  established  by  the 
hydraulic  principles  reported  in  1856  by  Henri  Darcy.  Darcy  discovered  that  the  flow  rate 
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through  porous  media  is  proportional  to  the  hydraulic  head  loss  and  inversely 
(Mroportional  to  the  length  of  the  flow  path.  Darcy's  law  can  be  stated  as 

Q  =  -KAVh=^qA  (2-5) 

where  Q  (mVsec)  is  the  volumetric  flow  rate.  A  (m^)  is  the  cross-sectional  area  of  the 
aquifer,  K  (m/sec)  is  the  hydraulic  conductivity,  Vh  (dimensionless)  is  the  gradient  in 
hydraulic  head,  and  q  (m/sec)  is  the  specific  discharge  or  flow  rate  per  unit  area.  The 
negative  sign  indicates  that  fluid  flows  in  the  direction  of  decreasing  hydraulic  head 
(Smith  and  Wheatcraft,  1993:6.7). 


Giound 


Figure  2.4  Hydraulic  Head  (Smith  and  Wheatcraft,  1993:6.7) 


Governing  Equations.  The  goal  of  groundwater  modeling  is  to  predict 
the  velocity  field  from  a  given  set  of  physical  conditions.  The  velocity  field  is  defined  by 
the  magnitude  and  direction  of  the  specific  discharge  rate.  The  development  of  the  flow 
equations  is  straightforward  because  all  of  the  flow  problems  are  developed  from  the 
same  fundamental  principle  of  conservation  of  fluid  or  mass.  The  derivation  is 
traditionally  done  by  referring  to  a  cube  of  porous  material  that  is  large  enough  to  be 
representative  of  the  properties  of  the  porous  medium  and  yet  is  small  enough  so  that  the 
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change  of  head  within  the  volume  is  relatively  snudl.  This  cube  of  pwous  material  is 
known  as  a  representative  elementary  volume  or  REV.  Its  volume  is  equal  to  AxAyAz 
(see  Figure  2.S). 

U 


AY 


Figure  Representative  Elementary  Volume  (Anderson  and  Woessner,  1992:16) 

The  general  conservation  equation  can  be  exf^essed  as 

m,  (2-6) 

where  m^  is  the  rate  of  change  in  mass  storage,  is  the  rate  of  mass  output,  and  is 
the  rate  of  mass  input.  Considering  the  x-direction  only,  the  change  in  mass  in  terms  of  q 
can  be  expressed  by 

=[q,(x+Ax)-q,(x)]AyAz  (2-7) 

where  q,(x-»-ZU)  is  the  x-conqwnent  of  discharge  at  (x-f  Ax),  q^(x)  is  the  x-component  of 
recharge  at  x,  and  AyAz  is  the  cross  sectional  area  of  the  REV  in  the  x-direction.  Similar 
reasoning  for  the  y-  and  z-directions  produces 

9,  ( Jf )]  AyAz + [q,{y + Ay)  -  q,  (y )]  AxAz + 

[q,(z  +  Az)“q,(z)]AxAy-RAxAyAz  (2-8) 

where  R  is  a  volumetric  source/sink  rate. 
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Using  the  (tefinition  of  specific  storage,  the  rate  of  change  in  sUMrage  in  the  REV  is  given 
by 

AV  Ah 

nt.-^^^S—AxAyAz  (2-9) 

At  At 

where  AV  is  the  change  in  fluid  volume,  is  the  specific  storage,  and  Ah  is  the  change  in 
hydraulic  head.  Combining  Equatimis  2-8  and  2-9,  dividing  through  by  AxAyAz,  and 
taking  the  limit  as  the  volume  goes  to  zero,  yields  the  final  form  of  the  water  balance 
equation: 


_c  » 

‘  dt  9jc  hy  9z 


(2-10) 


Because  the  specific  discharge  rate,  q,  cannot  be  measured  directly.  Equation  (2-10)  can 
not  be  used  to  compare  analytical  results  to  observed  field  data.  However,  hydraulic  head 
can  be  measured  directly  in  the  field  by  a  piezometer  and  the  gradient  can  be  estimated. 
Darcy's  law  can  then  be  used  to  determine  the  specific  discharge.  Darcy's  law  in  three 
dimensions  is 


_  „3h 
ft 


(2-11) 


Finally,  a  general  and  useful  form  of  the  governing  equation  for  hydraulic  head  can  be 
developed  by  substituting  Equation  (2-11)  into  Equation  (2-10)  yielding; 


(2-12) 


Every  problem  of  flow  that  can  be  modeled  is  described  by  one  or  more  equatitms  of  the 
form  found  in  Equation  (2-12)  (Anderson  and  Woessner,  1992:16-18). 
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Groundwater  Modeling 

Overview.  Hydrologists  are  often  called  upon  to  [nedict  the  response  of  an 
aquifer  system  to  proposed  changes.  The  concern  may  be  about  the  future  spatial 
distribution  of  water  levels,  water  quality,  or  cost  of  a  unit  volume  of  water  supplied  to  a 
customer.  The  tool  most  often  used  to  make  these  predictions  is  a  groundwater  flow 
model  of  the  aquifer  system. 

In  the  early  70's,  hydrologists  used  physical  models  such  as  sand  tanks  to  simulate 
groundwater  flow  directly.  Today,  most  hydrologists  use  mathematical  models  to 
simulate  groundwater  flow  indirectly  by  means  of  the  governing  equations.  A  conceptual 
model  of  the  aquifer  system  includes  descriptions  of  the  physical  processes,  hydraulic 
properties,  and  boundary  conditions  considered  relevant  to  the  particular  subsurface 
formation  under  study.  Once  formulated,  the  conceptual  model  is  then  translated  into  a 
mathematical  model.  Because  of  the  complexity  of  groundwater  systems,  most  models 
are  solved  numerically  rather  than  analytically.  Finite-difference  and  finite-element 
methods  are  the  two  most  commonly  used  techniques  for  solving  mathematical  models  in 
the  field  of  hydrology  (Anderson  and  Woessner,  1992:20).  The  main  features  of  these 
methods  are: 

1 .  The  solution  is  sought  for  the  numerical  values  of  state  variables  (hydraulic  heads) 
only  at  discrete  points  rather  than  their  continuous  variations. 

2.  The  partial  differential  equations  found  in  the  governing  equations  are  replaced  by  a 
set  of  algebraic  equations. 

3.  The  solution  is  obtained  for  a  specified  set  of  numerical  values  of  the  various  model 
parameters  rather  than  as  general  relationships. 

4.  Computer  codes  are  required  to  simultaneously  solve  the  very  large  number  of 
resulting  equations  (Bear  and  Verruijt,  1990:16). 
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Model  Calibration.  After  the  mathematical  model  has  been  developed  it  must 
be  calibrated.  Model  calibration  is  the  process  of  adjusting  tl^  model's  parameters  to 
obtain  a  reasonable  match  between  observed  site-specific  data  (calibration  targets)  and 
results  computed  by  the  model  (Walton,  1992:35).  In  effect,  the  model  is  calibrated  by 
reproducing  a  set  of  historical  data  with  some  acceptable  level  of  accuracy.  For 
groundwater  flow  models,  the  calibration  procedure  is  generally  carried  out  by  varying 
estimates  of  hydraulic  properties  from  a  set  of  initial  values  until  the  best  flt  of  calculated 
results  to  observed  data  is  achieved.  Model  calibration  is  synonymous  with  paran^ter 
estimation  and  is  often  referred  to  as  solving  the  inverse  problem. 

Model  calibration  is  one  of  the  most  difficult  and  time  consuming  aspects  of 
groundwater  modeling  (Water  Science  and  Technology  Board,  1990:225).  In  an  article 
on  model  calibration  techniques,  Barry  Power  cites  a  case  that  required  500  man-hours 
and  80  simulation  runs  to  estimate  one  parameter  in  finite-element  model  with  222  nodes 
(Power  and  Bams,  1993:9).  Some  of  the  reasons  parameter  estimation  is  so  difficult  are: 

1 .  Too  few  observations  are  available  to  compute  stable  estimates  of  statistics  such  as 
mean  and  variance. 

2.  Results  of  point  sampling  are  often  biased  broause  a  large  amount  of  data  does  not 
necessarily  allow  computation  of  nearly  true  or  effective  values  of  a  parameter  and  its 
variance.  For  example,  permeability  values  from  core  analyses  often  are  not 
representative  of  regional  values,  because  flow  through  large  fractures  is  not 
reproduced  by  core  analyses. 

3.  Observed  values  are  subject  to  numerous  sources  of  error  such  as  mismeasuring  water 
levels,  clogging  of  the  slots  or  screen  in  the  measuring  devices,  and  inaccurate 
reporting  (Cooley  and  Naff,  1989:4). 

The  techniques  for  parameter  estimation  are  numerous  and  range  from  simple 
graphical  curve  fitting  to  complex  statistical  estimation  algorithms.  Excellent  reviews  of 
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tlM  various  procedures  can  be  found  in  (Yeh,  1986)  and  (Carrera,  1988).  Histcaically,  the 
techniques  have  been  divided  into  two  major  categories:  (1)  manual  trail-and-error 
procedures  and  (2)  automated  or  optimization  methods  that  minimize  an  objective 
function. 

Manual  Trial>and*EiTor  Calibration.  In  manual  trial-and-error 
calibration,  the  parameters  are  assigned  initial  values  based  on  the  limited  information 
available.  During  calibration,  the  parameter  values  are  adjusted  sequentially  in  an  attempt 
to  match  simulated  heads  to  the  calibration  targets.  The  modeler  makes  the  parameter 
adjustments  based  on  his  or  her  expertise  and  knowledge  of  the  model  and  the  area  being 
simulated.  No  matter  how  the  method  is  applied,  it  has  several  inherent  deficiencies: 

1 .  It  typically  requires  tens  to  hundreds  of  model  runs  (Anderson  and  Woessner, 
1992:232). 

2.  No  methodology  exists  to  guarantee  that  simulations  will  proceed  in  a  direction  that 
could  lead  to  an  optimal  set  of  parameter  estimates. 

3.  It  is  difficult  to  determine  if  the  best  set  of  parameter  estimates  has  been  reached. 

4.  The  parameters  are  usually  calibrated  one  at  a  time,  overlooking  possible  interaction 
effects. 

5.  The  method  is  labor  intensive  (therefore  expensive),  frustrating  (therefore  often  left 
incomplete)  and  subjective  (therefore  biased)  (Carrera  and  Neuman,  1986:199). 

Manual  trail-and-error  was  the  first  technique  to  be  used  and  is  still  the  preferred  method 
in  practice  despite  its  many  drawbacks  (Keidser  and  Rosbjerg,  1991:2219). 

Automated  Calibration.  To  address  some  of  the  above  deficiencies, 
researchers  have  developed  several  automated  calibration  techniques.  Automated 
calibration  is  performed  using  specially  developed  codes  that  use  either  a  direct  or 
indirect  approach  to  solve  for  the  parameter  estimates  (Anderson  and  Woessner, 
1992:233).  In  a  direct  approach,  the  unknown  parameters  are  treated  as  dependent 
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variables  in  the  governing  differential  equations  and  heads  are  treated  as  independent 
variables.  To  use  this  approach,  the  head  values  must  be  specified  at  all  points  within  the 
system.  Since  heads  are  known  only  at  points  where  there  are  observation  wells,  the  head 
at  all  other  points  must  be  estimated.  For  this  reason  and  because  direct  solutions  are 
prone  to  instability',  the  direct  method  has  received  little  attention  and  is  not  often  used. 
The  main  advantage  of  the  direct  method  is  that  it  does  not  require  repeated  simulation 
runs. 

The  indirect  approach  is  similar  to  performing  manual  trial-and-error  calibrations  in 
that  the  simulation  is  executed  repeatedly.  However,  rather  than  adjusting  the  parameters 
based  on  the  modeler's  opinion,  some  objective  function  is  used  to  make  the  adjustments 
to  the  parameters.  The  objective  function  is  usually  some  measure  of  how  close  the 
calculated  heads  are  to  the  observed  heads.  Many  of  the  objective  functions,  also  called 
calibration  criteria,  can  be  written  in  a  general  form  as 

minimize  llh-hll  (2-13) 

A 

where  h  is  the  set  of  observed  heads  and  h  is  the  set  of  calculated  heads  (Lx)aiciga  and 
Church,  1990:645).  The  double  vertical  bars  represent  a  norm  or  measure  of  agreement. 
Some  of  the  calibration  criteria  that  have  been  used  in  the  past  are  (1)  mininuzation  of  the 
sum  of  squared  deviations  (i.e.,  least  squares  criterion),  (2)  minimization  of  the  maximum 
absolute  deviation,  and  (3)  minimization  of  the  sum  of  absolute  deviations.  Least  squares 
is  one  of  the  most  commonly  used  criteria  primarily  because  of  its  computational 
convenience.  Methods  used  to  optimize  (minimize)  the  objective  function  include  (1) 
linear  programming,  (2)  quadratic  programm  'ag,  (3)  Gauss-Newton  method,  and  (4) 
gradient  search  methods.  Most  studies  have  used  either  a  Gauss-Newton  or  gradient 
search  method  (Yeh,  1986:98). 

>  A  solution  is  instable  if  small  changes  (or  errors)  in  he'^.d  values  cause  large  changes  in 
the  parameter  estimates. 
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Automated  techniques  provide  a  systematic  approach  to  model  calibration.  They 
usually  require  fewer  simulation  runs,  are  less  subjective,  and  provide  statistical 
information  on  the  accuracy  of  the  parameter  estimates.  Despite  these  advantages,  they 
are  still  not  widely  used.  Perhaps  the  primary  reason  is  inertia  on  the  part  of  groundwater 
modelers.  However,  some  have  criticized  automated  techniques  for  failing  to  recognize, 
or  treat  adequately,  the  problems  of  nonidentifiability^,  nonuniqueness^,  and  instability 
(Carrera  and  Neuman,  1986:200). 

One  strategy  which  some  of  the  automated  techniques  touch  on  but  do  not  seem  to 
take  full  advantage  of  is  response  surface  methodology. 

Response  Surface  Methodology 

Response  surface  methodology  (RSM)  comprises  a  group  of  statistical  techniques  for 
empirical  model  building  and  model  exploitation.  RSM  (or  "The  methodology"),  first 
introduced  by  Box  and  Wilson  in  1951  and  later  developed  by  Box,  Hunter,  and  Draper, 
consists  of:  (1)  design  of  experiments,  (2)  regression  analysis,  and  (3)  model  exploitation 
or  optimization  (Cornell,  1990: 1).  One  of  the  primary  goals  of  RSM  is  to  find  the  best 
value  of  the  estimated  response.  In  the  case  of  parameter  estimation  for  groundwater 
models,  the  response  would  usually  be  the  residual  sum  of  squares  between  observed  and 
calculated  heads,  and  the  goal  would  be  to  find  the  set  of  parameter  values  which  produce 
the  minimum  error.  Figure  2.6  shows  an  example  of  a  sum  of  squares  error  surface  where 
the  model  has  two  parameters  (P,  and  Pj)- 


parameter  is  said  to  be  nonidentifiable  if  the  model  output  is  not  sensitive  to  it.  If  a 
parameter  is  nonidentifiable,  then  regardless  of  the  value  assigned  to  the  parameter,  the 
model  will  produce  the  same  output. 

^Even  if  all  the  parameters  are  identifiable,  the  solution  to  the  minimization  problem  may 
be  nonunique.  This  is  usually  associated  with  multiple  local  minima. 
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Figure  2.6  Sum  of  Squares  Error  Surface  (Box  and  others,  1978:463) 


Experimental  Design.  The  first  step  in  setting  up  an  experimental  design  is 
deciding  which  input  variables  (^,)  to  study.  Usually,  the  researcher  has  numerous 
possible  input  variables  to  consider.  In  applying  RSM  to  groundwater  flow  models,  the 
decision  involves  which  of  the  model  parameters  to  estimate  as  part  of  the  calibration 
process.  Often  for  convenience  and  computational  efficiency,  the  input  variables,  are 
transformed  into  coded  or  standardized  variables,  x,.  Generally,  the  following  linear 
transformation  is  used: 

(2-14) 

where  is  the  center  of  the  region  of  interest  and  5,  is  the  half  range  of  the  region.  The 
next  step  in  setting  up  a  design  determines  the  levels  of  the  input  variables  included  in  the 


design.  For  groundwater  flow  models,  a  reasonable  range  of  values  for  certain  parameters 
can  be  determined  from  observed  data.  For  example,  if  an  aquifer  is  mada  up  of  fine 
sands,  a  reasonable  range  for  the  levels  of  porosity  would  be  0.20-0.45  (see  Table  (2.1)). 
The  final  step  specifies  the  experimental  arrangement  or  design.  One  of  the  basic  designs 
is  a  two-level  factorial  (2*)  design.  In  such  a  design,  each  of  the  k  variables  occurs  at  just 
two  levels  (±1  in  coded  space).  Factorial  designs  in  general  have  the  following  useful 
properties: 

1 .  They  allow  numerous  comparisons  to  be  made  and  therefore  facilitate  model  creation 
and  criticism. 

2.  They  provide  estimates  whose  variances  are  as  small,  or  nearly  as  small,  as  those 
provided  by  any  other  design  occupying  the  same  space. 

3.  They  give  rise  to  simple  calculations  (Box  and  Draper,  1987: 106). 

Two-level  factorial  designs  are  often  used  to  estimate  models  with  main  effects  and  n- 
factor  interactions.  To  build  a  quadratic  model,  a  3*  factorial  arrangement  can  be  used, 
but  the  more  efficient  central  composite  designs  (CCDs)  and  Box-Behnken  designs  are 
generally  preferred.  A  central  composite  design  consists  of 

1 .  the  2*  vertices  of  a  ^-dimensional  cube  where  the  factor  levels  are  coded  so  that  the 
design  center  is  at  (0,0, . .  .,0).  The  values  of  the  coded  variables  in  this  factorial 
portion  are  (±14:1.  •  •  -.±1); 

2.  the  center  point,  no=(0,0, . .  .,0)  and; 

3.  the  2k  vertices  (±a,0,0, . .  .,0),  (04:0,0,  — 0), . . .,  (0,0, . .  .,04:a)  of  a  J^-dimensional 
star  (Cornell,  1990:52). 

Figure  2.7  shows  a  central  composite  design  for  three  variables. 

One  feature  often  incorporated  into  a  CCD  is  rotatability.  In  a  rotatable  design,  the 
accuracy  of  prediction  of  the  response  depends  only  on  its  distance  from  the  center  of  the 
design,  producing  a  variance  function  that  is  spherical  or  nearly  spherical.  Rotatability  is 
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achieved  by  assigning  impropriate  a-values.  A  general  fcmnula  used  to  achieve 
rotatability  in  a  full  facttMial  design  is: 


where  A:  is  the  number  of  factors,  is  the  number  of  observations  at  the  cube  points,  and 
is  the  number  of  observations  at  the  star  points  (Box  and  Draper,  1987:488). 
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Figure  2.7  Central  Composite  Design  (Cornell,  1990:53) 

Another  feature  often  incorporated  into  experimental  designs  is  orthogonality.  In  an 
orthogonal  two-level  design,  each  column  in  the  design  matrix  is  orthogonal  to  every 
other  column.  This  design  results  in  a  diagonal  variance-covariance  matrix  and  produces 
a  smaller  joint  confidence  region  than  a  nonorthogonal  design  (Box  and  Dither, 

1987:78).  Because  of  the  axial  points  in  a  central  composite  design,  the  column 
corresponding  to  the  constant  term  will  not  be  orthogonal  with  those  corresponding  to  the 
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quadratic  terms.  For  second  order  designs,  orthc^onality  refers  to  the  absence  of 
correlation  between  the  quadratic  terms. 

Box-Behnken  designs  are  a  subset  of  the  class  of  3^  designs.  By  leaving  out  the  comer 
points  from  the  fiiU  3*  design,  they  require  fewer  runs.  Box-Behnken  designs  are  often 
used  because  they  require  only  three  levels  for  each  variable  and  not  five  as  in  a  CCD 
(Cornell,  1990:60). 

Regression  Analysis.  After  selecting  a  design,  the  groundwater  model  can  be 
solved  and  the  residual  sum  of  squares  computed  for  each  design  point.  Regression 
analysis  techniques  can  then  be  applied  to  this  data  to  develop  an  enq>irical  model  of  the 
error  surface.  For  a  discussion  of  regression  analysis  techniques  including  the  method  of 
least  squares  and  analysis  of  variance,  see  any  standard  text  on  regression  analysis 
(Draper  and  Smith,  1981;  Montgomery  and  Peck,  1982). 

Optimizatioii.  The  third  area  of  RSM  is  searching  for  optimal  conditions.  To 
accomplish  this  task,  it  is  useful  to  know  what  type  of  surface  the  model  represents.  For 
second  order  models  with  only  one  variable,  a  simple  curve  is  produced.  Figure  2.8 
shows  a  curve  with  a  minimum.  If  the  model  had  included  two  variables,  some  of  the 
possible  surfaces  include:  a  simple  minimum  (or  maximum),  a  saddle,  a  stationary  ridge, 
or  a  rising  ridge  (see  Figure  2.9). 
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Figure  2.8  Simple  Curve  for  Second  Degree  Polynomial  in  One  Variable 


Figure  2.9  Ridge  Systems  for  Second  Degree  Polynomials  in  Two  Variables 

(Box  and  Draper,  1987:347) 
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For  second-order  models  with  three  variables,  only  the  iso-response  contours  associated 
with  the  fitted  surface  can  be  examined.  Figure  2.10  shows  some  of  these  contours.  It 
was  expected  the  iso-response  contours  of  the  residual  sum  of  squares  error  surface  would 
look  something  like  Figure  2.10a. 


Figure  2.10  Contour  Systems  for  Second  Degree  Polynomials  in  Three  Variables 

(Box  and  Draper,  1987:350) 
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There  are  two  basic  methods  used  to  search  for  optimal  conditicms:  (1)  a  calculus 
based  q>i»oach,  and  (2)  canonical  analysis.  In  the  calculus  based  iqq[>roach,  the  partial 
derivative  of  the  function  is  taken  with  respect  to  each  factor  and  each  partial  derivative  is 
set  equal  to  zero.  A  stationary  point  (if  one  exists)  can  then  be  found  by  solving  the 
resulting  set  of  equations.  The  nature  of  the  surface  at  the  stationary  point  can  be 
determined  by  examining  the  partial  second  derivatives.  The  partial  second  derivatives 
are  presented  in  the  Hessian  matrix: 


1 

¥ 

axjx* 

dxl 

If  x^Hx  ^  0  Vx  then  the  fimction  is  concave  and  will  have  a  simple  maximum  at  its 
stationary  point.  If  x'HElx  ^  0  Vx  then  the  function  is  convex  and  will  have  a  sinq)le 
minimum  at  its  stationary  point.  If  neither  of  the  above  is  true,  then  a  simple  maximum 
or  minimum  does  not  exist. 

The  other  method  used  to  analyze  second-order  models  is  canonical  analysis.  A 
canonical  analysis  involves  rotating  the  coordinate  axes  to  remove  all  cross-product  terms 
and  when  appropriate,  translating  the  coordinate  axes  to  coincide  with  the  stationary 
point.  Performing  a  canonical  analysis  simplifies  the  model,  allows  local  optimum  to  be 
easily  identified,  and  makes  the  description  of  the  surface  relatively  straightforward.  The 
details  of  canonical  analysis  are  available  in  several  references  (Box  and  draper,  1987; 
Khuri  and  Cornell,  1987;  Box  and  others,  1978) 
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m.  Results  and  Analysis 


Introduction 

The  primary  objective  of  diis  research  was  to  detennine  if  resprmse  surface  methodology 
(RSM)  could  be  used  as  a  parameter  estimatirm  technique  in  grou^^Klwater  flow  modeling. 
To  test  this  hypothesis,  an  atteiiq>t  was  made  to  calibrate  an  existing  two-dimrasioiud, 
steady-state  flow  model  using  the  RSM  approach.  The  model  was  reviewed  and  tested  on 
three  different  coiiq)uter  platforms.  Simulation  results  were  compared  to  the  output 
obtained  by  the.  developer  of  the  model.  The  available  field  data  was  examined  and  a 
calilM^on  criterimi  established.  Three  of  the  model’s  hydraulic  parameters  were  then 
selected  for  calibration  and  an  experimental  design  built.  Using  the  method  of  least 
squares,  a  regression  model  of  the  error  surface  was  developed  and  analyzed. 

LandfiU  Model 

Examining  the  efficacy  of  using  RSM  as  a  parameter  estimation  technique  for 
groundwater  flow  models  required  a  model  to  calibrate.  Rather  than  developing  a  new 
model,  this  research  conducted  experiments  on  an  existing  model. 

Model  Description.  The  particular  model  used  in  this  study  simulated 
groundwater  flow  and  solute  transport  for  a  cross-section  of  landfill  10  located  on 
Wright-Patterson  Air  Force  Base,  Ohio.  The  cross-section,  referred  to  as  JDb  in  the 
developer's  report,  is  oriented  in  a  generally  east-west  direction  and  extends 
approximately  1 100  feet  from  test  pit  WP-1JF10-TP12  to  monitoring  well  WP-LFIO- 
MW04  (see  Figure  3.1).  As  shown  by  the  lines  of  equal  potentiometric  head,  the 
groundwater  flow  within  this  area  was  nearly  aU  vertical.  The  porous  matrix  was 
composed  mostly  of  fine  sand,  gravelly  sand,  and  clay. 
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Figure  3.1  Cross<Section  JD  of  Landfill  10  (OSRI  Report,  1993:U*19) 

The  numerical  model  SUTRA  (Saturated'linsaturated  Transport)  was  used  to  simulate 
the  groundwater  flow  md  solute  transport  of  the  landfill.  SUTRA  is  a  computer  program 
developed  by  the  United  States  Geological  Survey  (USGS)  and  written  in  ANSI- 
STANDARD  FORTRAN-77.  The  version  of  SUTRA  used  in  this  study,  V-0690-2D, 
was  released  in  June  of  1990.  The  program  employs  a  two-dimensional  hybrid  finite- 
element  and  integrated  finite-difference  method  to  q)proximate  the  governing  differential 
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equations  (Voss,  1984:3).  This  approximation  method  divides  the  cross-sectional  aquifer 
domain  into  a  single  layer  of  contiguous  blocks,  called  finite  elements.  Flow  parameters 
and  variables  which  vary  from  point  to  point  in  the  aquifer  are  approximated 
elementwise,  nodewise,  or  cellwise  (see  Figure  3.2).  Elementwise  approximation  implies 
that  a  parameter  has  a  constant  value  over  each  element,  although  it  may  differ  from 
element  to  element.  Nodewise  approximation  results  in  a  continuous  surface  of  values 
with  linear  change  between  adjoining  nodes.  Cellwise  approximation  is  similar  to 
elementwise  approximation  but  each  cell  is  centered  on  a  node,  not  on  an  element.  The 
grid  for  cross-section  JDb  (Figure  3.3)  contained  1320  elements  and  1416  nodes. 
Constructing  a  finite  element  grid  is  a  time-consuming  process  and  one  of  the  reasons  an 
existing  model  was  used  in  this  study. 


Figure  3J2  Graphical  Representation  of  Elements,  Cells,  and  Nodes 
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Figure  3  J  Cross-section  JDb;  Finite  Element  Mesh 


Model  Execution.  SUTRA  requires  two  input  files  for  operation.  The  primary 
input  file  contains  all  the  data  necessary  for  the  simulation  except  the  initial  conditions. 
The  second  input  file  contains  the  initial  conditions  of  pressure  and  concentration.  As  the 
primary  output,  SUTRA  provides  fluid  pressures  and  solute  concentrations,  as  they  vary 
with  time,  everywhere  in  the  simulated  subsurface  system  (Voss,  1984:3).  The  model 
used  in  this  investigation  was  run  in  steady-state  mode  and  therefore  the  fluid  pressures 
did  not  vary  with  time.  The  model  also  calculated  solute  concentrations,  but  these  values 
were  not  relevant  to  this  study.  Copies  of  the  primary  input  file  and  output  file  for  the 


landfill  model  were  available  for  use  in  this  study.  The  initial  ccmditions  input  file  was 
not  available,  but  was  created  from  the  output  file  since  the  output  file  also  lists  the  initial 
conditions. 

SUTRA  can  be  compiled  and  executed  under  most  operating  systems  and  on  most 
computers.  The  developers  of  the  landfill  model  used  the  Lahey  fortran  compiler  and  a 
486  microcomputer.  The  Lahey  compiler  was  not  available  for  use  in  this  study.  The 
model  was  compiled  under  Microsoft  Fortran  version  S.  1  and  execution  was  attempted  on 
both  a  386  and  486  microcomputer.  On  both  machines,  the  program  would  not  execute 
properly.  The  model  was  successfully  compUed  and  executed  on  a  VAXA^S  machine 
and  under  the  UNIX  operating  system.  The  output  from  the  VAXA^S  and  UNIX  runs 
matched.  However,  there  were  significant  differences  between  the  results  obtained 
locally  and  the  output  supplied  by  the  developer  of  the  model.  The  pressure  differences  at 
the  1416  nodes  ranged  from  zero  to  8.9E03  N/m^.  Some  possible  explanations  for  these 
discrepancies  include  the  different  compilers  and  computer  platforms  used  and  the 
possibility  of  different  initial  conditions.  However,  these  differences  should  not  affect  the 
study  since  the  baseline  will  be  the  results  of  the  local  runs. 

Calibration  Criterion 

The  developer  of  the  model  used  manual  trial'and-error  and  graphical  matching 
techniques  to  calibrate  the  model.  The  developer  stated. 

Water-level  data  available  from  monitoring  wells  in  the  vicinity  of  Landfills  8  and  10 
were  used  to  aid  in  calibrating  the  numerical  models.  Hydraulic  parameters, 
particularly  the  degree  of  anisotropy  in  hydraulic  conductivity  values  and  the  rates  of 
recharge,  were  adjusted  in  the  model  until  water  levels  simulated  by  the  model 
resembled  water  levels  observed  in  the  monitoring  wells.  . . .  results  can  only  be 
tentative  until  rigorous  parameter  calibration  and  model  validation  is  performed. 
Herein,  the  model  is  used  to  simulate  plausible  hypotheses  and  aid  in  understanding 
the  physics  of  the  system.  Provided  that  the  uncertainty  of  model  parameters  is 
recognized,  this  understanding  can  be  helpful  in  assessing  the  risk  of  contaminant 
migration.  (OSRI  Report,  1993:5-49,U-7) 
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No  quantitative  information  on  the  closeness  of  observed  data  to  calculated  results  was 
provided. 

The  most  commonly  used  calibration  targets  are  observed  hydraulic  heads.  Head 
measurements  from  four  observation  wells  were  taken  on  three  separate  occasions  over 
an  eight  month  period  (see  Table  3.1).  The  three  observations  taken  at  each  well  were 
averaged  to  produce  a  single  calibration  target  per  well.  It  was  assumed  these  average 
values  represented  steady-state  conditions.  No  information  was  available  on  the  error 
associated  with  the  head  measurements. 


Table  3.1  Head  Observations 


Hydraulic  Head  (ft) 

Well  Number 

7  May  91 

19  Sen  91 

14  Jan  92 

Average 

Ol-GW-3 

904.05 

899.98 

906.40 

903.48 

WP-LF10-MW03A 

818.68 

818.53 

818.61 

818.61 

WP-LF10-MW04A 

797.80 

796.44 

796.58 

796.94 

WP-LF10-MW04B 

801.26 

799.66 

799.66 

800.19 

After  examining  the  model,  it  was  discovered  that  the  pressures  at  the  nodes 
corresponding  wells  WP-LF10-MW04A  and  WP-LF10-MW04B  were  set  to  constant 
values  as  boundary  conditions  for  the  model.  Consequently,  only  the  observations  for 
wells  Ol-GW-3  and  WP-LF10-MW03A  were  available  as  calibration  targets.  Since 
SUTRA  calculates  fluid  pressures  instead  of  heads,  the  fluid  pressures  for  the  two 
remaining  observations  were  determined  by 

p  =  (/»-z)pg  (3-1) 
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where  h  (m)  is  the  observed  hydraulic  head,  z  (m)  is  the  elevation  above  datum,  p  (kg/m^) 
is  the  fluid  density,  and  g  (m/sec^)  is  the  acceleration  due  to  gravity.  For  these 
calculations,  the  values  for  fluid  density  (1000  kg/m^)  and  gravitational  acceleration  (9.8 
m/sec^)  were  the  same  as  those  used  in  the  landfill  model.  The  elevation  values  used 
were  274.015  m  and  249.326  m  for  wells  Ol-GW-3  and  WP-LF10-MW03A  respectively. 
These  values  were  estimated  from  the  cross-sectional  diagram  shown  in  Figure  3.1.  The 
drawing  is  to  scale  and  the  estimates  were  obtained  by  using  the  middle  of  the  well 
symbol  as  a  reference  point.  The  resulting  approximations  seemed  reasonable  when 
compared  to  the  finite-element  mesh  which  gives  elevation  and  range  values  for  each 
node.  The  drawing  is  to  scale  and  The  results  of  converting  the  observed  hydraulic  heads 
to  fluid  pressures  appear  in  Table  3.2. 


Table  3.2  Observed  Fluid  Pressure 


Weil  Number 

Fluid  Pressure  (N/m^) 

Ol-GW-3 

1.33810563E+04 

WP-LF10-MW03A 

1.81336727E+03 

There  are  several  objective  functions  or  criteria  that  could  be  used  to  determine  the 
best  parameter  estimates.  Among  these  are  rrunimization  of  the  maximum  absolute 
deviation,  minimization  of  the  sum  of  absolute  deviations,  and  minimization  of  the  sum 
of  squares  of  the  deviations  (or  least  squares).  Because  of  its  widespread  use  and 
computational  convenience,  it  was  decided  to  use  the  least  squares  method.  The  equation 
for  calculating  the  residual  sum  of  squares  over  n  observation  points  is 

«SS  =  X(p,-P,)’  (3-2) 

i=I 
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where  p,  is  the  observed  pressure  (calibration  target)  and  p,  is  the  pressure  calculated  by 
the  model. 

To  establish  a  quantitative  value  for  how  close  the  model  was  initially  calibrated,  the 
developer's  output  was  compared  to  the  observed  fluid  pressures.  Table  3.3  shows  the 
comparison.  The  RSS  indicated  that  while  the  overall  results  may  match  the  observed 
water-level,  individual  comparisons  may  show  large  discrepancies. 


Table  33  Quantitative  Value  of  Calibration 


Well/Node 

Observed  Fluid 

Developer’s 

Residual  Sum  of 

Number 

Pressure  (N/m^) 

Results  (N/m2) 

Squares 

Ol-GW-3/244 

1.338 10563E+04 

-2.38676270E+04 

8.9116592E+09 

WP-LF10-MW03A/539 

1.81336727E+03 

8.85554844E+04 

Input  Variables 

The  model  contains  eleven  main  hydraulic  parameters  that  must  be  specified:  fluid 
density,  porosity,  permeability  (longitudinal  and  transverse),  fluid  viscosity,  porous 
matrix  compressibility,  fluid  compressibility,  residual  saturation,  parameters  of  the 
saturation-pressure  and  permeability  relationships,  and  rate  of  recharge  (OSRI  Report, 
1993:11-25,  U-29).  Some  of  these  parameter,  such  as  fluid  density,  fluid  viscosity, 
porous  matrix  compressibility,  fluid  compressibility,  and  residual  saturation,  can 
generally  be  adequately  determined  from  either  the  literature  or  from  field  and  laboratory 
experiments.  Values  for  the  others  must  usually  be  determined  through  a  calibration 
process.  To  demonstrate  the  RSM  technique,  three  parameters  (porosity,  transverse 
permeability,  and  rate  of  recharge)  were  selected  as  parameters  to  be  calibrated.  This 
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choice  seemed  reasonable  since  these  were  the  same  parameters  the  developer  of  the 
landfill  model  calibrated. 

Porosity.  The  porous  matrix  of  the  cross-sectional  layer  being  modeled  consists 
of  inorganic  clays,  gravelly  clays,  sandy  clays,  silty  clays,  lean  clays,  inorganic  silts,  very 
fine  sands,  rock  flour,  well-graded  sands,  and  gravelly  sands  (see  Figure  3.1).  As  shown 
in  Table  2.1,  typical  ranges  of  porosity  for  these  sediments  are:  clays  (0.40-0.60),  silts 
(0.35-0.50),  fine  sands  (0.20-0.45),  and  course  sands  (0.15-0.35).  The  developers  of  the 
model  used  0.40  as  the  value  for  porosity  (OSRI  Report,  1993:11-25).  This  study 
examined  porosity  values  fix>m  0.05  to  0.75. 

Transverse  Permeability.  Transverse  permeability,  when  different  from 
longitudinal  permeability,  gives  an  indication  of  the  degree  of  anisotropy  of  the  system. 
The  value(s)  for  transverse  permeability  cannot  be  determined  directly  from  field 
sampling.  The  developers  of  the  model  initially  assumed  an  isotropic  system  (1:1  ratio  of 
longitudinal  to  transverse  permeability),  but  later  calibrated  the  model  with  an  anisotropic 
ratio  of  100:1  (OSRI  Report,  1993:11-30).  This  study  examined  transverse  permeabilities 
that  produced  anisotropic  ratios  ranging  from  500: 1  to  55.6: 1. 

Rate  of  Recharge.  According  to  Figure  2.2  the  average  amount  of  yearly 
rainfall  reaching  the  groundwater  in  southwest  Ohio  is  four  inches  (not  including  run-off). 
When  the  developer  of  the  model  assumed  an  isotropic  system,  the  rate  of  recharge  was 
calibrated  to  an  unreasonable  value  of  18.6  inches  per  year  (OSRI  Report,  1993:11-30). 
When  the  model  v/as  run  as  an  anisotropic  system  ( 100:1),  the  rate  of  recharge  was 
calibrated  to  four  inches  per  year  matching  the  expected  value.  The  actual  parameter 
value  entered  into  the  model  is  computed  by  averaging  the  yearly  rate  over  the  surface 
area  (approximately  1 100  square  feet).  The  inflow  is  then  specified  at  each  of  the  surface 
nodes.  For  example,  the  four  inches  per  year  equates  to  2.0E-07  kg/sec  for  each  of  the 
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surface  nodes.  This  study  examined  rates  of  recharge  ranging  from  2.03E-09  kg/sec  to 
3.98E-07  kg/sec,  which  equate  to  0.04  and  7.96  inches  per  year,  respectively. 

Coded  Variables.  The  final  calibrated  value  for  each  parameter,  as  determined 
by  the  model  developer,  was  used  as  the  center  point  for  the  variables  in  the  design 
region.  The  three  variables  were  transformed  into  coded  space  by; 


^,-0.4  ^,-0.01  ^  ^j-2.00E-07 

■21 - ,  =  -22 - ,  and  or,  =  — - 

0.208  ^  0.00476  1.177E-07 


(3-3) 


Experimental  Design 

Since  a  sum  of  squared  residuals  was  used  as  the  calibration  criteria  for  this  study  (see 
Equation  (3-2)),  a  design  that  would  allow  the  creation  of  a  quadratic  model  was  needed. 
A  central  composite  design  (CCD)  was  chosen  over  a  31^  design  because  CCDs  are  more 
efficient,  requiring  fewer  simulation  runs.  Table  3.4  shows  the  CCD  design  in  both 
uncoded  and  coded  form.  The  design  is  a  composite  of  a  two-level  factorial  design 
represented  by  the  eight  cube  points  (±1,  ±1,  ±1),  combined  with  the  center  point  (0, 0, 0), 
and  augmented  with  six  star  or  axial  points  (±1.682, 0, 0),  (0,  ±1.682, 0),  (0, 0,  ±1.682). 
The  cube  points  allow  the  estimation  of  first  order  and  two-way  interaction  effects  while 
the  center  and  star  points  allow  the  estimation  of  pure  quadratic  effects.  The  star  points 
were  chosen  using  Equation  (2- IS)  to  make  the  design  rotatable.  Since  SUTRA  is  a 
deterministic  model,  only  one  run  was  conducted  at  each  of  the  design  points. 

Regression  Analysis 

The  results  of  executing  the  groundwater  flow  model  at  each  of  the  design  points  are 
shown  in  Table  3.5.  Table  3.6  lists  the  residual  sum  of  squares  and  the  sum  of  signed 
differences.  The  first  item  to  note  is  the  large  values  for  the  residuals  (on  the  order  of  10^ 
and  10'°).  Upon  obtaining  these  values,  the  calculations  involving  the  calibration  targets 
were  verified  and  found  to  be  correct. 
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Tables^  Central  Composite Dedgn 


Uncoded 

Coded 

Porosity, 

Transverse 
Permeability,  ^ 

Rate  of 
Recharge,  ^ 

Pwosity,  X, 

Transvase 
Permeability,  x, 

Rate  of 
Recharge,  x^ 

0.192 

0.00524 

8.23E-08 

-1 

-I 

-1 

0.608 

0.00524 

8.23E-08 

1 

-1 

0.192 

0.01476 

8.23E-08 

1 

-1 

0.608 

0.01476 

8.23E4)8 

1 

1 

-1 

Cube 

0.192 

0.00524 

3.18E-07 

-1 

1 

Points 

0.608 

0.00524 

3.18E-07 

1 

1 

0.192 

0.01476 

3.18E-07 

1 

1 

0.608 

0.01476 

3.18E-07 

1 

1 

1 

0.400 

0.01000 

2.00E-07 

0 

0 

0 

Center  Point 

0.050 

0.01000 

2.00E-07 

-1.682 

0 

0 

0.750 

0.01000 

2.0(®-07 

1.682 

0 

0 

0.400 

0.00200 

2.00E-07 

0 

-1.682 

0 

Axial 

0.400 

0.01800 

2.00E-07 

0 

1.682 

0 

Points 

0.400 

0.01000 

2.03E-09 

0 

0 

-1.682 

0.400 

0.01000 

3.98E-07 

0 

0 

1.682 

Tables^  Results  from  CCD 


Run 

Coded  Variables 

Fluid  Pressure 

Porosity, 

XI 

Transverse 
Permeability,  x-y 

Rate  of 
Recharge,  x-t 

Node  244 

Node  539 

1 

-1 

-2.24948402E+04 

7.11774274E+04 

2 

1 

-1 

-2.24948402E-i4)4 

7.11774274E+04 

3 

-1 

1 

-1 

-2.45766239E-f04 

9.96608946E404 

4 

1 

1 

-2.45766239E+04 

9.96608946E404 

5 

-1 

1 

-2.19215516E-f04 

7.11849601E404 

6 

1 

-1 

1 

-2.19215516E+04 

7.11849601E+04 

7 

1 

1 

-2.42890995E-fO4 

9.96709079E+04 

8 

1 

1 

1 

-2.42890995E404 

9.96709079E+04 

9 

0 

0 

0 

-2.36630370E+04 

8.84435243E404 

10 

-1.682 

0 

0 

-2.36630370E+04 

8.84435243E-K)4 

11 

1.682 

0 

0 

-2.36630370E-t04 

8.84435243E-f04 

12 

0 

-1.682 

0 

•2.48210152E404 

1.05434369E+05 

13 

0 

1.682 

0 

•1.88335580E+04 

5.12087129E-fO4 

14 

0 

0 

-1.682 

-2.335345 17E+04 

8.845 13403E404 

15 

0 

0 

1.682 

-2.3986401 8E-f04 

8.8435 1885E+04 
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Table  3.6  RSS  and  Sum  of  Signed  Deviatioiis 


Run 

Residual  Sum  of 
Squares  (RSS) 

Sum  of  Signed 
Deviations 

1 

6.0984S278SE+09 

-3.348816369E+04 

2 

6.0984S2785E+09 

-3.348816369B»04 

3 

1.1014924088E+10 

-5.988984718E+04 

4 

1.1014924088E+10 

-5.988984718E+04 

5 

6.05869201 3E409 

-3.406898499E+04 

6 

6.058692013E409 

-3.406898499E404 

7 

1.0995138893E+10 

-6.018738489E+04 

8 

1.0995138893E+10 

-6.018738489E4O4 

9 

8.877048953E+09 

-4.958606379E-fO4 

10 

8.877048953E+09 

-4.958606379E404 

11 

8.877048953E+09 

-4.958606379E+04 

12 

1.2196710264E+10 

-6.541893029E404 

13 

3.477681542E409 

•  1.71 8073 138E-f04 

14 

8.855562446E-f09 

•4.990346509E^ 

15 

8.899666835E409 

-4.925436319E+04 

Some  possible  explanations  for  the  large  discrepancies  between  observed  and  calculated 
pressures  include: 

1 .  The  boundary  conditions  of  the  model  may  not  be  correctly  specified. 

2.  The  two  nodes  from  the  finite-element  mesh  (244  and  539)  may  not  correspond 
exactly  to  the  locations  of  the  observation  weUs. 

The  second  item  that  can  be  easily  seen  from  Table  3.5  is  that  porosity  has  no  effect  on 
the  RSS.  Comparing  the  runs  where  the  porosity  value  goes  fiom  low  (-1)  to  high  (-i-l) 
while  the  other  parameter  remain  unchanged  clearly  shows  that  porosity  has  no  effect  on 
the  output  of  the  model.  This  may  be  a  case  of  nonidentifiablity  as  discussed  in  chq)ter  2. 
Despite  the  unusually  large  residuals  and  porosity  not  effecting  the  output,  an  en^)irical 
model  of  the  error  surface  was  developed.  As  discussed  earlier,  it  was  believed  a 
quadratic  model  of  the  form 

gCjf.P)  —  Po  ■*’Pl'*l  "^^2^2  ■^p3'*3  ■*'Pll^I  '^P22'*2  ■*'P33'*3 

PljXjXj +P13X1X3 +P23X2X,  (3-4) 
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would  be  needed  to  adequately  represent  the  error  surface.  The  statistical  software 
package,  SAS,  was  used  to  analyze  the  RSS  data.  The  regression  i»ocedure  was  run  with 
the  input  variables  in  their  coded  form.  The  best  model  produced  was 

y  =  8.892616142  10’  +  2.1516519152  10’xj-3.65644160  10*;tj"  (3-5) 

As  expected,  the  porosity  variable  was  not  included  in  the  model.  Also,  the  rate  of 
recharge  was  determined  to  have  no  significant  effect  on  the  output,  leaving  only  the 
transverse  permeability  variable  in  the  model.  The  associated  analysis  of  variance  for  the 
model  is  shown  in  Table  3.7.  Even  though  this  is  the  "best"  model  for  the  data,  as 
evidenced  by  the  large  F  and  R-Square  values,  the  mean  square  error  and  the  standard 
errors  for  the  estimates  were  extremely  large.  Also  the  coefficient  of  variation  (C.V.), 
which  expresses  the  ratio  of  the  root  mean  square  error  to  the  mean  response,  was  high. 
All  these  indicators  imply  that  the  model  has  a  high  degree  of  variability  and  is  of  low 
precision. 


Table  3.7  Analysis  of  Variance  for  Equation  3*5 


Sum  Of 

Mean 

Source 

DF 

Squares 

Square 

F  Value 

Prob>F 

Model 

2 

8.8042932E19 

4.4021466E19 

8575.668 

0.0001 

Error 

12 

6.1599585E16 

5.1332987E15 

C  Total 

14 

8.8104532E19 

RootMSE  7.1646704E07  R-Square  0.9993 

Dep  Mean  8.5596789E09  Adj  R-Square  0.9992 

C.V.  0.83703 


Variable 

DF 

Parameter 

Estimate 

Standard 

Error 

T  for  HO: 
Param=0 

Prob>m 

Intercept 

1 

8.892616142E09 

2.66463 18E07 

333.728 

0.0001 

Xj 

1 

2.516519152E09 

1.9386560E07 

129.807 

0.0001 

X* 

i2 _ 

1 

-3.65644160E08 

2.1062304E07 

-17.360 

0.0001 
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A  gn4>h  of  the  regression  noodel  rei»esented  by  Equation  (3-5)  is  shown  in  Figure  3.4. 
The  area  between  the  two  vertical  lines  is  the  experimental  region.  This  gnq>h  shows  that 
the  optimal  value  for  transverse  permeability  was  not  in  the  design  region  and  that 
lowering  the  value  of  transverse  permeability  fvoduces  lower  values  for  the  residual  sum 
of  squares.  The  graph  also  indicates  that  the  anisotropic  ratio  of  100:1,  which  the 
developer  of  the  landfill  model  considered  to  be  close  to  optimal,  was  not  warranted. 


Figure  3.4  Plot  of  Regression  Model  (Equation  3-5) 

Given  the  poor  utility  of  the  results,  the  basic  assunqitions  of  this  study  were  reviewed. 
The  first  assumption  was  that  the  calibrated  parameter  values  determined  by  the 
developer  of  the  landfill  model  were  close  to  optimal.  However,  Table  3.3  showed  that. 
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at  least  for  the  two  observation  points  available,  a  significant  difference  existed.  Another 
assumption  was  that  porosity  and  rate  of  recharge  would  effect  the  model's  response.  The 
results  in  Table  3.5  and  the  absence  of  these  parameters  from  the  regression  model 
demonstrated  that  they  had  little  or  no  effect.  Based  on  the  review  of  these  assuiiq)tions, 
the  lack  of  additional  data,  and  the  low  i»ecision  of  the  regression  model,  it  was  decided 
further  investigation  with  the  RSM  techniqiM  would  be  pointless  in  this  case.  While  no 
definite  conclusions  about  the  effectiveness  of  RSM  as  a  parameter  estimation  technique 
could  be  made  from  the  results  of  this  study,  the  potential  still  remains.  The  next  step  in 
this  process  would  be  to  test  the  technique  on  a  standardized  or  hypothetical  model  and 
possibly  using  a  screening  design  to  determine  which  parameters  have  the  most  effect  of 
the  model's  response. 
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rV.  Conclusioiis  and  Recommendations 


Conclusions 

The  goal  of  this  research  was  to  detennine  if  response  surface  methodology  (RSM)  could 
be  used  to  estimate  groundwater  flow  parameters.  This  study  was  motivated  by  the  need 
for  a  systematic  {q)proach  to  groundwater  flow  model  calibration.  Currently,  manual 
trial-and-error  techniques  are  the  most  commonly  used  calibration  methods. 

An  attempt  was  made  to  calibrate  an  existing  two-dimensional,  steady-state  flow 
model  using  RSM.  The  developer  of  the  landfill  model  used  a  gnqihical  technique  to 
calibrate  several  of  the  model's  hydraulic  parameters.  The  developer's  approach  consisted 
of  matching  predicted  water-table  levels  to  observed  water-table  levels.  The  developer 
stated  the  levels  were  close  but  no  quantitative  information  on  the  closeness  of  fit  was 
provided.  The  RSM  approach  was  applied  using  the  developer’s  calibrated  parameter 
values  as  the  starting  point  of  the  investigation  and  a  residual  sum  of  squares  was  used  as 
the  calibration  criterion.  Following  the  developer's  lead,  porosity,  transverse 
permeability,  and  rate  of  recharge,  were  the  three  hydraulic  parameters  which  were  to  be 
estimated.  After  executing  the  landfill  model  at  each  of  the  design  points,  the  residual 
sum  of  squares  were  calculated.  The  residuals  had  an  order  of  magnitude  of  10^  and 
indicating  a  significant  difference  between  the  observed  and  calculated  fluid  pressures. 
Also  the  regression  model  of  the  error  surface  revealed  a  large  amount  of  variability  in  the 
data.  Only  the  transverse  permeability  parameter  was  included  in  the  regression  model 
indicating  that,  at  least  in  the  area  of  the  design  region,  porosity  and  rate  of  recharge  had 
little  or  no  effect  on  the  landfill  model's  output. 

This  investigation  started  with  the  assumption  that  the  parameter  values  determined  by 
the  landfill  model  developer  were  close  to  optimal  and  that  the  three  hydraulic  parameters 
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selected  for  calibration  would  be  significant.  The  assumptions  turned  out  not  to  be  valid. 
Because  of  these  bad  assumptions,  low  precision  of  the  regression  model,  and  lack  of 
additional  data,  the  investigation  was  terminated.  No  definite  conclusions  about  the 
effectiveness  of  RSM  as  a  parameter  estimation  technique  could  be  made  from  the 
results. 

Recommendations 

The  following  are  some  recommendation  for  extending  this  research: 

1 .  Rather  than  using  a  groundwater  flow  model  of  a  real  system,  the  RSM  technique 
could  be  tested  against  a  hypothetical  or  synthetic  model.  An  example  of  using  a 
hypothetical  model  can  be  found  in  an  article  by  Chu,  Strecker,  and  Lettenmaier  (Chu 
and  others,  1987). 

2.  A  preliminary  screening  design  could  be  used  to  determine  which  hydraulic 
parameters  have  the  most  effect  on  the  response  surface  and  then  those  parameters 
used  for  further  study. 

3.  An  investigation  could  be  conducted  to  compare  the  results  of  using  different 
calibration  criteria,  such  as  minimization  of  the  absolute  differences,  minimization  of 
the  maximum  difference,  and  a  weighted  least  squares. 
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